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Abstract 

We propose a solution procedure for the finite element discretiza- 
tion of a static nonlinear compressible Mooney-Rivlin hyperelastic 
solid. We consider the case in which the boundary condition is a large 
prescribed deformation, so that mesh tangling becomes an obstacle for 
straightforward algorithms. Our solution procedure involves a largely 
geometric procedure to untangle the mesh: solution of a sequence of 
linear systems to obtain initial guesses for interior nodal positions for 
which no element is inverted. After the mesh is untangled, we take 
Newton iterations to converge to a mechanical equilibrium. The New- 
ton iterations are safeguarded by a line search similar to one used in 
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optimization. Our computational results indicate that the algorithm 
is much faster than a straightforward Newton continuation procedure 
and is also more robust (i.e., able to tolerate much larger deforma- 
tions) . 

1 The problem under consideration 

We consider the problem of solving for the deformed shape of a hyperelastic 
solid body under static loads. The continuum mechanical model under con- 
sideration has the following description [0]. Let Bq C R'^ be an undeformed 
solid body whose boundary is OBq. Here d, the space dimension, is 2 or 3. 
Assume boundary conditions (either displacement or traction, i.e., Dirichlet 
or Neumann) are given as follows. The boundary OBq is partitioned into 
two subsets Yd and r^r. A function (po : ^ R*^ specifies new-position 
(Dirichlet) boundary conditions. A second function to : Fat — R'^ specifies 
traction (Neumann) boundary conditions. Everything in this paper extends 
to the more general case that some coordinate entries are Neumann while 
others are Dirichlet at certain boundary points, but we limit the discussion 
to the special case that each boundary point is Dirichlet or Neumann in all 
d coordinates in order to simplify notation. Finally, the model requires a 
specification of the model's body forces, that is, a function h : Bq ^ R*^ that 
specifies the force of gravity and other forces on the body. 

The problem is to find a function (p : Bq —>■ that specifies the new 
position of the body. Let B denote (j){Bo). For a point X e Bq, let x = 
0(X). Let F be the deformation gradient, i.e., F = d(f)/dK. = rfx/rfX. It is 
assumed that -F(X) has a positive determinant for all X. The Green strain 
is defined to be = {F'^F — I)/2. Let scalar function \l/(-F) be the strain 
energy function, which is assumed to be a property of the material. For 
this paper, we assume that \E' depends only on two matrix invariants of E, 
namely J = det(F) = ^det{2E + I) and h = trace(F^F) = trace(2E + J). 
Further specializing the model, the strain energy is then taken to have the 
following form suggested by Ciarlet and Geymonat in P] for compressible 
Mooney-Rivlin materials 

^(F) = ^( - 1) - Q + /i) In J + ^(/i - 3) (1) 

where A, /i > are material parameters. For the d = 2 case, we assume that 
Bq is 3D but that the z-displacement is identically and that the x- and 
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^-displacements depend only on x and y; these are called plane strain as- 
sumptions. Thus, E has a last row and column of all zeros, and the Mooney- 
Rivlin formula in is applied to this E to come up with the strain energy 
function for the d = 2 case. The condition for static equilibrium (written in 
minimization form) is that 

[ ^{F(X))dV- I ph-(P{X)dV- I to-0(X)rfA (2) 

is minimized among all choices of (p that satisfy the Dirichlet boundary con- 
dition, i.e., that satisfy (/>(X) = 4>o(X.) for all X G F/j. 

This condition can be rewritten in variational form: for all admissible 
variations 6u, that is, functions in the space [H^{Bq)Y that vanish on Td, 

f ^ : Giad Su dV - [ ph-SudV-f to-SudA = 0, (3) 

J Bo J Bo Jtn 

where A : B = trace(^i?^) is used to denote the inner product of second- 
order tensors A and B. This model also applies to the case of linear elasticity 
with two changes in definitions. First, E = {{F - I)^ + {F - I))/2 in the 
case of linear elasticity. Second, \1/ = j E{i, j)"^ + | -^(^; 0)^; which 
can be written in terms of the two matrix invariants of E. 

We should mention that our method does not appear to depend so much 
on the specific details of the Mooney-Rivlin model, except for the In J term, 
which is quite important for our analysis. Since dv = JdV, where dv is the 
volume element of B and dV is the volume element of Bq, this logarithmic 
term resists infinite compression of the material: if a small positive volume 
of material in Bq shrinks to a 0- volume set in B, then this term causes the 
strain energy at those points to become infinite. 

We next describe the discretization of the problem under consideration 
pp. We assume that Bq is discretized with a mesh of triangles or tetrahedra. 
We assume that the discretization of 0, or alternatively the discretization 
of the displacement u = 0(X) — X, is piecewise linear, with the pieces of 
linearity being the mesh cells. (In Section |Sl we discuss extension of our 
method to piecewise quadratic displacements.) This assumption implies that 
u is determined by dm real numbers, namely, the values of u at nodes. Recall 
that d is the space dimension, and let m denote the number of non-Dirichlet 
nodes of the mesh. The finite element method finds the displacement u such 
that holds for all test functions 6u in the test function space. Here, the 
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test function space is the set of (5u's that are piecewise hnear and continuous 
and vanish on F^). The integral in is evaluated with a quadrature rule; 
we have used a 6-point formula having degree 4 precision from (7j for our 
quadrature in 2D and a 15-point formula having degree 5 precision from 
jHj for 3D. These quadrature rules were chosen because they have positive 
weights; the Gauss points are interior; and they possess full triangular or 
tetrahedral symmetry, respectively. It suffices to solve for the dm choices 
of 6u that compose the standard basis for the test function space. This yields 
a system of dm nonlinear equations for dm unknowns. 

The algorithmic question under consideration is how to robustly solve 
these nonlinear equations. In the next section, we give a summary of the 
mesh tangling issue and of our proposal to overcome it. The individual steps 
of our algorithm are then described in more detail in Sections |3] and |S| Our 
computational experiments are presented in Sections El and [7| Concluding 
remarks are presented in Section |H1 

2 Mesh tangling 

The standard method for solving a system of nonlinear equations is Newton 
iteration. It is well-known, however, that if the initial guess is far from the 
true solution, then Newton iteration will often diverge. 

In the case of hyperelasticity with large deformation, there is a specific 
obstacle that may cause divergence, namely, mesh tangling. The definition 
of this term is that a mesh is tangled if the value of J defined in the previous 
section is or negative in Bq. In the case of linear displacements, J is 
piecewise constant, and hence this condition can be verified with a finite 
number of determinant computations. The matter of checking for tangling in 
the piecewise quadratic case is more complicated and is discussed in Section|H] 
A solution with a tangled mesh is physically invalid. Indeed, the strain 
energy function is undefined in this case because of the presence of the term 
— (I + /i) In J. Note that although the strain energy function is undefined 
when J is negative, the Galerkin form (j21) is still well defined, which is an 
anomaly that we return to below. We assume that the given problem instance 
has a valid solution, i.e., there is a piecewise linear function u satisfying the 
boundary conditions, as well as Q for all test functions Su plus the condition 
that J > on every element. 

Even with this assumption, Newton's method will still often run into 
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problems because the mesh will become tangled on intermediate steps. For 
example, the starting point for Newton's method is often taken to be u = 
on every interior node. If the deformation of the boundary is large, then 
this starting point will have tangling among most of the elements that are 
adjacent to the boundary. 

To understand a difficulty posed by a tangled mesh, suppose that the 
strain energy has a single term \1'( J) = — (| + /x) In J (for | + /i > 0) on a 
single element and that there are no boundary constraints. If we treat J as 
the independent scalar variable, then Newton's method for minimizing this 
scalar function is 

j(*+i) = _ 
which simplifies in this case to 

j(.+i) = 2JW. 

For positive J^^\ this iteration produces a sequence of J's tending to +00. 
This is to be expected since the minimum of — (| + /i) In J is indeed at +cx3. 
On the other hand, for a negative J^'^\ this iteration tends to —00, which is 
physically invalid. 

The preceding analysis, although naive, seems to point to the following 
conclusion: Newton's method on the Galerkin form, when applied to a tan- 
gled mesh, has a natural tendency to make the tangling worse. We suspect 
that this fact is probably already known to experts in the field, although we 
have not been able to find it in the previous literature. 

Given the conclusion in the previous paragraph, it seems of paramount 
importance to avoid tangling. When Newton's method fails in computational 
mechanics, it is standard practice to try Newton continuation, that is, to 
apply the load in incremental steps and use the converged solution for one 
step as the Newton starting point for the next step. Continuation is described 
in more detail in Section El Continuation, however, addresses the tangling 
issue only in an indirect fashion and therefore is likely to be very inefficient. 
Our computational experiments confirm the inefficiency of continuation. 

We propose a new algorithm for getting around the mesh tangling ob- 
stacle. The basic idea is to first untangle the mesh using a much simplified 
mechanical model. Once the mesh is untangled, the true mechanical model is 
solved. "Untangling the mesh" means finding a (j) that satisfies the Dirichlet 
boundary condition and also satisfies J > 0. Our new algorithm, which we 
call UBN (for "untangling before Newton" ) consists of three steps. 
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1. First, we attempt to untangle the mesh with the iterative-stiffening 
algorithm, described in Section EJ Iterative stiffening builds on the 
FEMWARP algorithm from our previous work [T^. If the iterative 
stiffening algorithm cannot untangle the mesh, then UBN reports fail- 
ure to solve the problem. 

2. Else if iterative stiffening succeeds, then we take Newton iterations to 
solve Q. The starting point for Newton is an untangled mesh produced 
by step 1. No continuation is used. On the other hand, Newton's 
method is safeguarded using a line search described in Sectional which 
prevents the introduction of new tangling. The line search is based on 
a technique common in the interior point literature. 

3 Newton continuation 

In the case that direct use of Newton's method to find fails to converge, 
the standard alternative is Newton continuation, also known as applying the 
load in steps. This section briefly describes Newton continuation before we 
return to a description of UBN. Newton continuation is also used to make 
graphs of load-displacement relationships, but in this paper we assume that 
the problem under consideration is to determine a single final configuration 
rather than an entire loading path. 

The basic form of Newton continuation is quite straightforward: a se- 
quence of parameters = Ao<Ai<---<A7v = lis chosen, and a sequence 
of displacement vectors uq, ui, . . . , ujv is computed, in which for each k, 
is the solution to the discretized in the case that (0o(X) — X, b, to) are 
replaced by ■ (0o(X) — X, b, to). Solution uq (corresponding to absence 
of loads) is identically 0. Solution is found via Newton's method, where 
Ufc_i is used as the initial guess. The final deformed configuration is given 
by u^r since Aat = 1. Note also that it is possible to accept a low- accuracy 
(not fully converged) solution for when k < N since it is presumably not 
necessary to achieve high accuracy for intermediate results that are not part 
of the ultimate answer. 

In some cases, a straight linear parametrization of the load path (as in 
the previous paragraph) is not feasible. In this case, one must construct a 
nonlinear parametrization ((/)ArLp(X; A), bArip(X; A), to,ArLp(X; A)) with the 
property that 0lp(X; 0) = X while 0lp(X; 1) = 0o(X) and similarly for the 
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other load terms. Examples of nonlinear parametrizations are given later in 
the paper. 

The only remaining issue is how to select the sequence of A^'s. We use an 
adaptive rule defined as follows. Assume that there are no body forces and 
that the traction boundary conditions are all zero (i.e., "traction-free" sur- 
faces) . This means that the only loading term is the Dirichlet boundary con- 
dition. We form the deformed mesh Mfe_i after applying the displacements 
given by Uk-i to non-Dirichlet nodes and Dirichlet boundary conditions 
scaled by Xk-i, (i.e., the deformed position is given by X + Afc_i(0o(X) — X) 
in the case of linear parametrization) to Dirichlet nodes. Next, we compute 
a value of such that, if the boundary nodes in M^-i are further deformed 
to positions given by X -|- Afc(0o(X) — X), then no tetrahedron altitude will 
decrease by more than a factor of r], where rj < 1 is a. tuning parameter of the 
continuation algorithm. (In particular, the step is sufficiently small that the 
mesh will not tangle after the new boundary condition given by A^ is applied 
to M.) This adaptive strategy appears to work reasonably well, although we 
did encounter some robustness problems discussed in Section [7| 

The description in the previous paragraph assumed the special case of 
traction-free Neumann boundaries and absence of body forces. It is more 
difficult to use this adaptive technique when there are nonzero body forces or 
tractions since it is not obvious how to step these loads in a way that prevents 
tangling on each step. Therefore, our test cases involve only traction-free, 
body-force-free loading conditions. Since the focus of the paper is the UBN 
method, it represents a strengthening of our contention that UBN is usually 
better than the competing algorithm (continuation) since we limit our testing 
only to the case that seems well suited for continuation. 

4 Iterative stiffening for mesh untangling 

In this section, we describe our procedure called iterative stiffening for un- 
tangling a mesh. We take the original mechanical problem given by Q, 
and using the same boundary conditions and loads, we solve the equations 
of isotropic linear elasticity using piecewise linear (constant-strain) finite el- 
ements j]p. Note that these equations have the same material parameters 
(the Lame constants A and /x) as the Mooney-Rivlin model. Linear elastic- 
ity requires one linear system solve. If the deformed mesh (i.e., the mesh 
that arises from moving the vertices to their displaced positions) is untan- 
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gled, then iterative stiffening is finished. If not, then our iterative stiffening 
procedure locates all elements that are inverted in the deformed mesh and 
increases their stiffness by 50%. The linear elasticity model is now solved 
again. This procedure is repeated indefinitely until the mesh is untangled or 
an excessive number of iterations has passed. 

We have not found this precise version of iterative stiffening appearing in 
the previous literature, but it is related to ideas already in the literature. It 
could be regarded as an extension of FEMWARP, a finite element based mesh 
warping approach developed by the authors within the linear weighted lapla- 
cian smoothing (LWLS) framework JHI- One difference is that FEMWARP 
does not easily encompass the mechanical concept of traction boundary con- 
ditions. It is also related to a mesh warping method used for ALE solvers 
and described in Chapter 7 of [T]. 

In our preliminary version of the UBN method J7], the untangling was 
done using Opt-MS 4j rather than iterative stiffening. Opt-MS is an untan- 
gling algorithm that iteratively repositions interior nodes one at a time until 
the mesh is untangled. It solves a small linear-programming problem for each 
node to find the position for it that maximizes the minimum determinant of 
its neighbors. We found recently that iterative stiffening is more effective 
for use in UBN than Opt-MS. One possible reason is that it is difficult to 
implement traction boundary conditions in a natural way in Opt-MS. 

Note that iterative stiffening can be made particularly efficient by using 
matrix-updating. In particular, it is well-known (see, e.g. 0) how to update 
a Cholesky factorization of a symmetric positive definite matrix A after A 
has undergone a low-rank update. If the iterative stiffening procedure stiffens 
only a few elements per iteration (our test runs confirm that indeed there 
are usually only a few updates per step), then this can be implemented as a 
a low-rank update, which is potentially much more efficient than solving a 
new stiffness matrix from scratch. We did not implement matrix- up dating 
because the work for iterative stiffening was usually dominated by the solver 
part of the algorithm anyway. 

5 Newton Line Search 

Newton's method is often employed for solving nonlinear systems of con- 
tinuously different iable equations 0. Let / : R*^*" R'^™', continuously 
differentiable, be given. The task at hand is to find a u G R'^™ such that 
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/(u) = 0. Let uq G R'^"* be given. Then, at each iteration k, Newton's 
method solves 

J(ufe)sfc = -/(ufc), (4) 

where J denotes the Jacobian of /, for the Newton step, s^., and performs 
the following update 

Ufc+i = Ufc + Sfc. (5) 

If it becomes necessary to satisfy one or more additional inequality con- 
straints, it is possible to safeguard the Newton step with the introduction of 
a line search. Let ak denote the line search parameter. Then ak is chosen to 
be as large as possible such that < ajt < 1 and u^+i = + a^Sfc satisfies 
the constraint. 

It is often difficult to compute the value of that minimizes /(xfe + a^Sfc) 
and satisfies the constraint because / is often a highly nonlinear function. 
In addition, the optimal value of ak often produces steplengths that are too 
short in practice. Thus, it is common practice in interior point methods to 
derive heuristics for computing that allow for both ease of computation 
and larger steplengths [12]. One such heuristic is to choose so as to stay 
a fixed percentage away from the boundary. We employ this heuristic in our 
line search below. 

As was pointed out in Section |2l the mesh is tangled unless J = det(F) > 
0. Thus, we introduce a line search that enforces that J > on each iteration 
of Newton's method. In particular, we begin with J > on the zeroth itera- 
tion and choose the line search parameter such that J(ufe_|_i) > O.lJ(ufc) 
on each element so as to stay a fixed percentage away from the boundary for 
reasons discussed above. 

The following pseudocode algorithm shows how the line search parameter 
is determined. Let denote the number of elements in the mesh. Given 
a displacement vector u, it is straightforward for each element i = 1, . . . ,N 
to compute the deformation gradient F and its determinant J determined 
by this displacement on element i; we denote the resulting determinant by 
J{u,i). 

Let Uq be the value of the displacement (at non-Dirichlet nodes) returned 
by the previous iteration of our Newton/line search algorithm. Initially, the 
value of Uq is the output of the iterative stiffening algorithm. It is assumed 
that the mesh determined by the Dirichlet boundary conditions and by uq on 
non-Dirichlet nodes is untangled. Let s denote the Newton step determined 
from Uq via (j3)). 
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a = 1; 

hTi = l: N 
while true 

if J(uo + as, i) > J(uo) / 10 
break 

end 

a = a ■ 0.9; 

end 

end 

6 2D Experiments 

We designed a series of numerical experiments in order to test the robustness 
of UBN and to compare it to the standard Newton continuation algorithm. 
For all of the numerical experiments in this paper, we set the parameters in 
(P) as follows: A = (^^_^_J^^_2i,) and /i = 2(1+^0' ^^^^ E = 1 and z/ = 0.3. 

The termination criteria for the Newton loop in UBN and for the final 
step of Newton continuation was that ||-F||2 < 10~^°||-Fo||2, where Fq is the 
initial value (i.e., the value when all interior displacements are set to 0) of 
the load vector. For the Newton continuation steps prior to the final step, 
the termination criteria was that ||-F||2 < 10~^||-FfcJ|2; where F^.- is the initial 
value of the load vector at the beginning of major iteration k. 

The algorithms were implemented in Matlab. The current version of the 
solver does not allow specification of nonzero traction boundary conditions 
nor specification of body forces. There is no difficulty in principle in ex- 
tending the Matlab code for UBN to handle these types of loads, but new 
quadrature formulas would have to be written. As mentioned earlier, there 
are some conceptual difficulties in extending our adaptive method for select- 
ing Afc's for continuation to these types of loads. 

The linear solution operation in Matlab is quite highly optimized and is 
expected to compete well with a custom-written C or C++ linear solver. On 
the other hand, the matrix assembly process involves several nested Matlab 
loops and is therefore expected to be much slower than a C or C++ version. 
For this reason, wall-clock times derived from the Matlab code are not useful 
predictors of computational demands that would be observed with a C or 
C++ code. 

Instead, we measure the running time in terms of assembly/linear solve 
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steps. An assembly/linear solve (ALS) step consists of one stiffness matrix 
and load vector assembly operation followed by one sparse linear system 
solve. The Newton continuation method involves a sequence of Newton solve 
procedures, and each Newton solve is further subdivided into several ALS 
steps. The UBN method involves iterative stiffening iterations followed by 
a safeguarded Newton method. We count each iteration of iterative stiffen- 
ing as an ALS step. The assembly portion of the iterative stiffening ALS 
operation is not exactly the same as the assembly portion of Newton, since 
the former involves linear elasticity assembly whereas the latter involves non- 
linear tangent stiffness assembly. We ran both assembly codes on an older 
Windows machine running Matlab 5.3, which has a "flops" function built 
in that measures floating point operations. (Newer versions of Matlab lack 
this function.) From this experiment we determined that the number of op- 
erations for the two kinds of assembly are fairly close. Furthermore, both 
assembly operations are much less costly than the linear system solve. Note 
that the iterations of iterative stiffening would be considerably cheaper than 
an ALS step had we implemented low-rank corrections described in Section |3] 

The solver portion of UBN involves additional operations connected with 
the line search. We determined (again by running test cases in Matlab 5.3) 
that the line search requires a tiny number of operations in comparison to 
the solution of the linear equations. 

Thus, it is sensible to compare the running time of UBN to continuation 
by considering the total number of ALS steps required for either. 

In this section we describe our 2-dimensional test case, which is an annu- 
lar domain. The mesh was generated with Shewchuk's Triangle J3] and is 
illustrated in Fig. ^ It contains 182 nodes and 286 triangles. 

The boundary conditions used in this test case involve a rotation of the ex- 
terior boundary circle by / radians combined with moving the inner boundary 
by a factor / closer to the outer boundary (where / = means no motion and 
/ = 1 means that the inner boundary would coincide with the outer bound- 
ary). Values of / tried were 0.1, 0.3, 0.6, and 0.7. The resulting deformed 
meshes are illustrated in Fig. |21 

The number of ALS steps to compute these deformed meshes is given 
in Table H The columns of this table are as follows. The first column / is 
the amount of boundary deformation as described in the previous paragraph. 
The second column UBN-IS is the number of iterations of iterative stiffening 
required by UBN. The third column UBN-NM is the number of iterations 
of Newton's method required by UBN. The fourth column UBN-ALS is the 
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Figure 1: Annulus used for testing in this section. 



number of ALS steps required by UBN (and hence is the sum of the second 
and third columns). 

The remaining columns of the table report on results from the continu- 
ation algorithm. The fifth column is the number of major iterations (i.e., 
updates to the continuation parameter A) required by the continuation algo- 
rithm using the adaptive rule discussed in Section IHl with parameter r] = 1/3. 
The sixth column is the number of ALS steps required by continuation. The 
seventh and eighth columns are the same quantities when rj = 1.2. Note that 
?7 = 1.2 is a quite aggressive choice of stepsize for continuation, since any 
value of ?7 > 1 means that updating the boundary could cause an inversion. 
For this 2D test case, an aggressive choice of t] did not seem to hinder con- 
vergence, but the results for a large value of i] in 3D described in the next 
section are less favorable. 

For / = 0.8, neither method was able to find a solution. UBN's iterative 
stiffening did not untangle the mesh after the maximum number of iterations 
(400) had been reached, and continuation stalled at A = 0.93. 

It should be noted that a highly deformed mesh like the solution when / = 
0.7 is probably not physically valid because the finite element discretization 
is no longer an accurate approximation to the underlying PDE. Nonetheless, 
we include extreme cases like this because it is interesting to compare the 
two algorithms in limiting cases. The test results show that UBN is much 
faster than continuation for both modest and extreme deformations. 

Note that for continuation, the outer boundary motion is parametrized 
by 6, the rotation angle. Linear parametrization with respect to (x, y) co- 
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Figure 2: Deformed annulus meshes with / = 0.1, / = 0.3, / = 0.6, and 
/ = 0.7 respectively. 

ordinates would work poorly in this case because a linear deformation from 
the initial position of the outer boundary to the final position would cause 
the outer boundary to shrink in radius and then expand. 

Comparing the UBN-ALS and Contin-ALS columns of this table indi- 
cates that UBN is approximately 15 to 30 times more efficient than continu- 
ation when f] = 1/3. Continuation is 3 to 4 times faster when used with the 
larger value of rj but is still significantly slower than UBN. Other annulus 
deformation tests not reported here confirm that UBN is always far more 
efficient than continuation. 

We also wished to check whether the line search procedure built into UBN 
was ever active in order to determine whether it is an essential part of UBN. 
We found that it was active on about 30%-50% of the iterations for the larger 
values of deformation. 
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Table 1: Results of comparison of UBN to continuation for the 2D annular 
domain. See the text for meanings of the columns. 
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7 3D Experiments 

Our experiments in 3D consisted of two tetrahedral meshes called "Hook" 
and "Foam5," which were provided to us by P. Knupp [9 . The sizes of the 
meshes are as follows: Hook contains 1190 vertices and 4675 tetrahedra, and 
Foam5 contains 1337 nodes and 4847 tetrahedra. Hook is contained in a 
bounding box of size 54 x 40 x 95, while Foam5 is contained in a bounding 
box of size 11.3 x 5.5 x 6.6. 

In both cases, we applied Dirichlet boundary conditions to two of the 
boundary surfaces, leaving the rest traction-free. In both cases, the Dirichlet 
conditions are identically zero on one boundary surface and are pulling in a 
uniform direction on the other surface by three magnitudes. For Hook, the 
displacement sizes were 10, 20, and 40, whereas for foam they were 0.5, 2, 
and 5. Thus, we see that the applied displacements are on the same order as 
the size of the object, and therefore large deformations will result. Figure El 
shows the deformed and undeformed configurations of Hook for the maximum 
deformation of 5, while Fig. |3]shows the corresponding illustration of Foam5. 

The results of our tests of UBN versus Newton continuation on Hook and 
Foam5 are given in Table|21 As in the previous section, the unit for measuring 
running time is ALS steps. For these tests, straight linear parametrization 
was used for continuation. As before, our results on Hook indicate that UBN 
is 15 to 50 times faster than continuation when rj = 1/3, and 5 to 15 times 
faster for rj = 1.2. 

The continuation algorithm terminates when the increment in A became 
smaller than the prespecified minimum (0.0005); this happened in two of the 
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(c) id) 



Figure 3: The top row (a) and (6) shows the undeformed Hook mesh from 
two different viewpoints. The dark blue surface is traction free, while the 
Dirichlet boundary conditions are applied to the green and pink surfaces. 
The bottom row (c) and (d) show Hook after the maximum deformation of 
40 is applied. 
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(c) id) 

Figure 4: The top row (a) and (b) shows the undeformed FoamS mesh from 

two different viewpoints. The dark blue surface is traction free, while the 
Dirichlet boundary conditions are applied to the green surface and pink edge. 
The bottom row (c) and (d) show FoamS after the maximum deformation of 
5 is applied. 
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Table 2: Results of comparison of UBN to continuation for 3D domains. See 
the text for meanings of the columns. 



mesh 


applied displ. 




UBN 






Continuation 














V = 


1/3 


V = 


1.2 






IS 


NM 


ALS 


Majit 


ALS 


MajIt 


ALS 


Hook 


10.0 


1 


5 


6 


53 


107 


15 


32 


Hook 


20.0 


1 


6 


7 


105 


212 


30 


61 


Hook 


40.0 


1 


8 


9 


210 


421 


59 


119 


FoamS 


0.5 


1 


4 


5 


28 


36 


*** 


*** 


FoamS 


2.0 


1 


5 


6 






*** 


*** 


FoamS 


5.0 


1 


7 


8 






*** 


*** 



tests with FoamS when = 1/3 (indicated by ' — ' in the table). Apparently 
this is due to an extremely flat tetrahedron which, although it is does not 
become inverted, causes the heuristic used for adaptively incrementing A to 
take very conservative steps. 

The continuation algorithm also terminates when inverted elements re- 
main after one major iteration. This happened in three of the tests with 
Foam5 when rj = 1.2 as indicated by '***' in the table. This shows that, as 
expected, the aggressive choice of r] may be more prone to inverting elements. 

8 Conclusions 

In summary, we developed a robust solution method for solving nonlinear 
elasticity equations for hj^erelastic solids with large boundary deformations. 
The basic idea is to first untangle the mesh using purely geometric methods 
and second solve the mechanical model; thus, the algorithm was named UBN 
(for "untangling before Newton"). The first step of our algorithm is to at- 
tempt to untangle the mesh with iterative stiffening. Assuming the mesh is 
untangled, UBN takes safeguarded Newton steps to solve 0. 

We tested the robustness of UBN and compared it to the standard New- 
ton continuation algorithm. We demonstrated that UBN is much faster than 
the Newton continuation algorithm. It could be argued that continuation 
would be more competitive with UBN if only we had used a different strat- 
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egy for incrementing A^. This may be true, but it seems to us that there is no 
good universal fast method for choosing the A^. Our experiments indicate, 
for example, that a more aggressive algorithm for updating A is more prone to 
terminating early due to inverted elements. Even selecting the continuation 
path seems to be nontrivial (e.g., for the 2D annulus example, it was nec- 
essary to parametrize the Dirichlet boundary condition in polar rather than 
rectangular coordinates). In contrast, the UBN method does not require any 
such analogous problem- dependent decisions, and the only parameters of the 
algorithm involve termination criteria. 

As described so far, our method applies to finite elements in which the 
displacement field is piecewise linear over tetrahedra, but UBN could be ex- 
tended to piecewise quadratic displacements. The challenge with piecewise 
quadratic displacements is that checking for tangling is much more compli- 
cated, as J is not constant on the element. In particular, it is a function 
of both the displacement and the location on the element, which makes it 
difficult to determine when J > analytically. There are some separate nec- 
essary and sufficient conditions for element inversion in the literature. Let G 
be the Jacobian of F. Then one such necessary condition is that det((j) has 
the same sign (strictly positive or strictly negative) at some finite list of test 
points [Hj . In this case, we would test for inversion at the Gauss points used 
for numerical quadrature; however, it is still possible that folding could occur 
at the corners. A more complicated sufficient condition for invertibility in- 
volving the Bernstein-Bezier form of a polynomial is given in jT^. Checking 
that the sufficient condition is met requires running a linear programming 
algorithm. Salem, Canann, and Saigal have proposed sufficient conditions 
for quadratic triangles and tetrahedra in [TU], [TT], [T2], and |13]. 

Finally, it should be noted that there are some classes of problems that 
are not amenable to UBN. For example, consider a long cylinder in which one 
end is held at zero displacement, the other is rotated by 2ti radians, the long 
side-surface is traction-free, and there are no body forces. For this problem, 
UBN would compute a solution where all displacements are identically zero. 

The difficulty with this problem instance lies in the fact that for this ge- 
ometry and specification of boundary conditions, there are an infinite number 
of continuum solutions (namely, rotation of one endcap by any integer mul- 
tiple of 27r). To distinguish one solution from another requires additional 
information beyond boundary conditions. It is not clear what form that ad- 
ditional information ought to take. (For this particular problem instance, 
the additional data is the particular sought-after multiple of 27r, but it is not 
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clear what specification is needed in general.) It is also not clear how UBN 
should be extended in order to use such additional information. 

9 Acknowledgements 

The authors wish to thank Patrick Knupp of Sandia National Laboratories 
for providing us with the 3D test meshes. They benefited from many helpful 
conversations with Katerina Papoulia of Cornell University. 

References 

[1] T. Belytschko, W. Liu, and B. Moran. Nonlinear Finite Elements for 
Continua and Structures. Wiley, 2000. 

[2] P. Ciarlet and G. Geymonat. Sur les lois de comportement en elasticite 
non-lineaire compressible. C. R. Acad. Sci. Paris Ser. II, 295:423-426, 
1982. 

[3] J. E. Dennis, Jr. and R. B. Schnabel. Numerical Methods for Uncon- 
strained Optimization and Nonlinear Equations. Society for Industrial 
and Applied Mathematics, 1996. 

[4] L. Freitag and P. Plassmann. Local optimization-based simplicial mesh 
untangling and improvement. Int. J. Numer. Meth. Eng., 49:109-125, 
2000. 

[5] G. H. Golub and C. F. Van Loan. Matrix Computations. John Hopkins 
University Press, 3rd edition, 1996. 

[6] G. Holzapfel. Nonlinear Solid Mechanics: A Continuum Approach for 
Engineering. Wiley, 2000. 

[7] T. J. R. Hughes. The Finite Element Method: Linear Static and Dy- 
namic Finite Element Analysis. Dover Publishers, 2000. 

[8] P. Keast. Moderate-degree tetrahedral quadrature formulas. Comput. 
Method. Appl. M., 55:339-48, 1986. 

[9] P. Knupp. Personal Communication, October 2002. 



19 



[10] A. Salem, S. Canann, and S. Saigal. Robust quality metric for quadratic 
triangular 2D finite elements. Trends in Unstructured Mesh Generation, 
220:73-80, 1997. 

[11] A. Salem, S. Canann, and S. Saigal. Mid-node admissible spaces for 
quadratic triangular 2D arbitrarily curved finite elements. Int. J. Nu- 
mer. Meth. Eng., 50:253-272, 2001. 

[12] A. Salem, S. Canann, and S. Saigal. Mid-node admissible spaces for 
quadratic triangular 2D finite elements with one edge curved. Int. J. 
Numer. Meth. Eng., 50:181-197, 2001. 

[13] A. Salem, S. Saigal, and S. Canann. Mid-node admissible space for 3D 
quadratic tetrahedral finite elements. Eng. Comput., 17:39-54, 2001. 

[14] M. Shephard, S. Dey, and M. Georges. Automatic meshing of curved 

three-dimensional domains: Curving finite elements and curvature- 
based mesh control. In I. Babuska, J. Flaherty, J. Hopcroft, W. Hen- 
shaw, J. Oliger, and T. Tezduyar, editors. Modeling, Mesh Genera- 
tion, and Adaptive Numerical Methods for Partial Differential Equa- 
tions. Springer Verlag, 1995. 

[15] J. Shewchuk. Triangle: Engineering a 2D quality mesh generator and 
Delaunay triangulator. In Proceedings of the First Workshop on Ap- 
plied Computational Geometry, pages 124-133, New York, NY, 1996. 
Association for Computing Machinery. 

[16] S. Shontz and S. Vavasis. An algorithm based on finite element weights 
for warping tetrahedral meshes. Archived by arxiv.org, cs.NA/0410045, 
2004. 

[17] S. M. Shontz. Numerical Methods for Problems with Moving Meshes. 
PhD thesis, Cornell University, 2005. 

[18] S. Vavasis. A Bernstein-Bezier sufficient condition for invertibility of 
polynomial mappings. Archived by arxiv.org, cs.NA/0308021, 2003. 

[19] S. J. Wright. Primal-Dual Interior- Point Methods. Society for Industrial 
and Applied Mathematics, 1997. 



20 



